20 October 2017

Why?

Why use R as a GIS?

1- Workflow

2- Quite efficient

3- Interface

Workflow

R for everything (almost)

  • import your data
  • format your data
  • analyze your data
  • visualize your data
  • export your data
  • create your own function/package

Efficiency

  • well-defined classes
  • many formats (+ convert different file)
  • topology operations
  • statistical analyses
  • visualization

Interface

Overview

Turning R into a powerful GIS

1- Define classes

  • sp
  • raster

2- Topology operations

  • rgeos

3- Import / export

  • rgdal

Turning R into a powerful GIS

5- Analyses

  • spdep
  • dismo

4- Creating and editing maps

  • graphics
  • ggplot
  • interactive plot

6- But

  • Geo-referencing …

Package sp

Classes

Classes / Functions Contents
Points list of points (set of coordinates)
SpatialPoints list of points + CRS
SpatialPointsDataPoints list of points + CRS + attribute table
Line a line (set of coordinates)
Lines list of lines
SpatialLines list of lines + CRS
SpatialLinesDataFrame list of lines + CRS + attribute table

Example: SpatialPointsDataPoints

long lat var1 var2
-80.69108 43.82453 0.3886067 4.807244
-81.89949 43.33342 0.6131389 1.598417
-80.04768 43.55549 0.8443004 8.640650
-81.21052 43.89272 0.9910437 4.709529
-80.72660 42.69218 -0.6817615 2.370249
-80.50302 42.86226 -0.8261933 7.134282

Example: SpatialPointsDataFrame

library(sp)
## Loading required package: methods
mysp <- SpatialPointsDataFrame(
  coords = mydata[,1:2],
  data = mydata[,3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Example: SpatialPointsDataFrame

class(mysp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
isS4(mysp)
## [1] TRUE
slotNames(mysp)
## [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

Example: SpatialPointsDataFrame

mysp@proj4string
## CRS arguments:
##  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
head(mysp@data)
##         var1     var2
## 1  0.3886067 4.807244
## 2  0.6131389 1.598417
## 3  0.8443004 8.640649
## 4  0.9910437 4.709529
## 5 -0.6817615 2.370249
## 6 -0.8261933 7.134281

Change projection: spTransform

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))
##            coordinates       var1     var2
## 1  (-8982490, 5408730)  0.3886067 4.807244
## 2  (-9117009, 5333529)  0.6131389 1.598417
## 3  (-8910867, 5367457)  0.8443004 8.640649
## 4  (-9040314, 5419220)  0.9910437 4.709529
## 5  (-8986444, 5236254) -0.6817615 2.370249
## 6  (-8961555, 5261955) -0.8261933 7.134281
## 7  (-9034590, 5390266) -1.6096683 5.141021
## 8  (-8934524, 5133405) -1.3334443 3.679047
## 9  (-9007348, 5265309)  1.0359543 1.289148
## 10 (-9125375, 5430049) -0.7524389 2.482117
## 11 (-8956996, 5260018) -0.6503393 5.434667
## 12 (-9052092, 5270655) -0.2185950 8.072032
## 13 (-9052073, 5338709) -0.7221735 4.921236
## 14 (-9072344, 5215846)  0.2134285 5.715557
## 15 (-8944403, 5162613) -0.3092620 4.673095
## 16 (-9032815, 5196734) -1.3373535 9.352931
## 17 (-8976914, 5232732)  0.2677742 5.277918
## 18 (-9058252, 5284466) -1.3750872 2.602788
## 19 (-8916892, 5373576) -0.4299518 4.955816
## 20 (-9026967, 5410482) -1.0852490 1.382881

Package raster

Classes

1- SpatialGrid (package sp)

2- SpatialGrid (package sp)

3- RasterLayer (package raster)

4- RasterStack (package raster)

5- RasterBrick (package raster)

Example:raster

library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
## 
##     extract
ras1 <- raster(matrix(runif(100*100,0,10),ncol=100,nrow=100),
    crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
    xmn=-82, xmx=+80, ymn=+42, ymx=+44)

Example:raster

ras1
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1.62, 0.02  (x, y)
## extent      : -82, 80, 42, 44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0006633392, 9.998927  (min, max)

Retrieving free GIS data: getData

head(getData("ISO3"))
##   ISO3        NAME
## 1  ABW       Aruba
## 2  AFG Afghanistan
## 3  AGO      Angola
## 4  AIA    Anguilla
## 5  ALA       Åland
## 6  ALB     Albania

Example:getData

## Country level:
mapBEL0 <- getData(name="GADM", country="BEL", path="./assets", level=0)
mapBEL1 <- getData(name="GADM", country="BEL", path="./assets", level=1)
tminW <- getData(name="worldclim", var="tmin", res=10, path="./assets")
mapBEL0
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 2.555356, 6.40787, 49.49722, 51.50382  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 68
## names       : OBJECTID, ID_0, ISO, NAME_ENGLISH, NAME_ISO, NAME_FAO,              NAME_LOCAL, NAME_OBSOLETE,            NAME_VARIANTS, NAME_NONLATIN, NAME_FRENCH, NAME_SPANISH, NAME_RUSSIAN, NAME_ARABIC, NAME_CHINESE, ... 
## min values  :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ... 
## max values  :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ...

Package rgdal

Drivers

library(rgdal)
## rgdal: version: 1.2-13, (SVN revision 686)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5
  1. writeOGR / writeGDAL: to write spatial objects
  2. readOGR/readGDAL: to read spatial files

Drivers

head(ogrDrivers())
##         name                     long_name write  copy isVector
## 1 AeronavFAA                   Aeronav FAA FALSE FALSE     TRUE
## 2 AmigoCloud                    AmigoCloud  TRUE FALSE     TRUE
## 3     ARCGEN             Arc/Info Generate FALSE FALSE     TRUE
## 4     AVCBin      Arc/Info Binary Coverage FALSE FALSE     TRUE
## 5     AVCE00 Arc/Info E00 (ASCII) Coverage FALSE FALSE     TRUE
## 6        BNA                     Atlas BNA  TRUE FALSE     TRUE

Export

writeOGR(mysp, dsn="./assets", layer="mypoints",
    driver="ESRI Shapefile", overwrite_layer=TRUE)

Import

mysp2 <- readOGR(dsn="assets/", layer="mypoints")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "mypoints"
## with 20 features
## It has 2 fields
## Roads find on http://www.diva-gis.org/Data
canroads <- readOGR(dsn="assets/", layer="CAN_roads")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "CAN_roads"
## with 16915 features
## It has 5 fields

Note: there are also functions to import/export raster in raster

Package rgeos

Load the package

library(rgeos)
## rgeos version: 0.3-25, (SVN revision 555)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE

Load the package

buf <- gBuffer(mapBEL0, width=0.5)
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
diff <- gDifference(gBuffer(mapBEL0, width=0.5), gBuffer(mapBEL0, width=0.1))
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
## Warning in gBuffer(mapBEL0, width = 0.1): Spatial object is not projected;
## GEOS expects planar coordinates

Package mapview

Import package

library(mapview)
## Loading required package: leaflet

NB: it uses the leaflet package.

Quick examples

mapview(mysp, cex = 'var2')@map

Quick examples

mapview(mapBEL1)@map

More

Editing a map

A very basic map – Shapefile

plot(mapBEL0)

A very basic map – Raster

class(tminW)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

A very basic map – Raster

plot(tminW)

Customize a map – Shapefile

plot(mapBEL0, border='grey15', col='#E6E6E6', lwd=1.6)
plot(mapBEL1, lty=2, lwd=0.9, add=T)
points(4.3513, 50.8471, pch=19, col="#27df9d")
text(4.3513, 50.8471, text="Brussel", pos=3)

Customize a map – Shapefile

Customize a map – Shapefile of roads

Part 2 – Showcases

Data frame with Ontario's lakes

  • the .csv file below is available in the repository
lakedf <- read.csv('assets/lakeOnt.csv')
head(lakedf)
##   X LAKEID LONGITUDE_DD LATITUDE_DD Area_ha
## 1 1 L32173    -80.03361    49.90417   174.9
## 2 2 L21001    -91.68278    48.68028   130.2
## 3 3 L44001    -84.19667    46.91861   279.6
## 4 4 L25002    -90.43111    48.16639   114.1
## 5 5 L44002    -84.30528    47.06000   137.3
## 6 6 L21003    -91.34528    48.24806  2993.9

Make it a SpatialPointDataFrame

library(sp)
lakesp <- SpatialPointsDataFrame(
  coords = lakedf[,3:4],
  data = lakedf[,c(1,4)],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Your SpatialPointDataFrame

lakesp
## class       : SpatialPointsDataFrame 
## features    : 546 
## extent      : -95.06111, -76.01194, 43.05722, 54.52056  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## variables   : 2
## names       :   X, LATITUDE_DD 
## min values  :   1,    43.05722 
## max values  : 546,    54.52056

Retrieve altitude for Canada

altCAN <- getData(name="alt", country="CAN", path="./assets/")
altCAN
## class       : RasterLayer 
## dimensions  : 4992, 10620, 53015040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/kevcaz/Github/Tutorials/mapsWithR/docs/assets/CAN_msk_alt.grd 
## names       : CAN_msk_alt 
## values      : -108, 5800  (min, max)

Retrieve altitude for Canada

Retrieve provinces' boundaries

bouCAN <- getData(country='CAN', level=1, path="./assets/")
bouCAN
## class       : SpatialPolygonsDataFrame 
## features    : 15 
## extent      : -141.0069, -52.61889, 41.67693, 83.11042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,     TYPE_1, ENGTYPE_1, NL_NAME_1,                                       VARNAME_1 
## min values  :        1,   42, CAN, Canada,    1, Alberta,  CA.AB,    NA,    10,   Province,  Province,          ,                                                 
## max values  :       15,   42, CAN, Canada,   13,   Yukon,  CA.YT,    NA,    62, Territoire, Territory,          , Yukon Territory|Territoire du Yukon|Yukon|Yuk¢n

Retrieve provinces' boundaries

bouCAN@data[,1:8]
##    OBJECTID ID_0 ISO NAME_0 ID_1                    NAME_1 HASC_1 CCN_1
## 1         1   42 CAN Canada    1                   Alberta  CA.AB    NA
## 2         2   42 CAN Canada    2          British Columbia  CA.BC    NA
## 3         3   42 CAN Canada    3                  Manitoba  CA.MB    NA
## 4         4   42 CAN Canada    4             New Brunswick  CA.NB    NA
## 5         5   42 CAN Canada    5 Newfoundland and Labrador  CA.NF    NA
## 6         6   42 CAN Canada    6     Northwest Territories  CA.NT    NA
## 7         7   42 CAN Canada    7               Nova Scotia  CA.NS    NA
## 8         8   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 9         9   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 10       10   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 11       11   42 CAN Canada    9                   Ontario  CA.ON    NA
## 12       12   42 CAN Canada   10      Prince Edward Island  CA.PE    NA
## 13       13   42 CAN Canada   11                    Québec  CA.QC    NA
## 14       14   42 CAN Canada   12              Saskatchewan  CA.SK    NA
## 15       15   42 CAN Canada   13                     Yukon  CA.YT    NA

Canada – Ontario

bouCAN[11,]
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -95.15598, -74.34605, 41.67693, 56.86972  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,   TYPE_1, ENGTYPE_1, NL_NAME_1,    VARNAME_1 
## min values  :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada 
## max values  :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada

Canada – Ontario

Canada – Ontario

plot(crop(altCAN, bouCAN[11,]@bbox))
par(mar=c(0,0,0,0))

Canada – Ontario

Canada – Ontario elevation

altONT <- rasterize(bouCAN[11,], crop(altCAN, bouCAN[11,]@bbox), mask=TRUE)
mapview(altONT)+mapview(lakesp)

Canada – Ontario elevation

Canada – Ontario elevation

par(mar=c(0,0,0,0))
plot(altONT)
plot(lakesp, add=T)

Canada – Ontario elevation

Resources

Useful links